package pl.softech.demo2d.physics.fluid;

import pl.softech.demo2d.physics.fluid.grid.Cell;

public class NSDiscretization {

	private double dx, dy;
	private double dt;
	private double dx2;
	private double dy2;
	private double dx2dy2;
	private double dxdy;
	private double gx, gy;

	/**
	 * Lepkość kinematyczna
	 * 
	 * http://en.wikipedia.org/wiki/Viscosity#Kinematic_viscosity:_.CE.BD_.3D_.
	 * CE.BC_.2F_.CF.81
	 * 
	 */
	private final double nu;

	public NSDiscretization(double dx, double dy, double dt, double gx, double gy, double nu) {
		setDx(dx).setDy(dy).setG(gx, gy);
		this.dt = dt;
		this.nu = nu;
	}

	public double getDx() {
		return dx;
	}

	public double getDy() {
		return dy;
	}

	public double getDt() {
		return dt;
	}

	public void setDt(double dt) {
		this.dt = dt;
	}

	public double getGx() {
		return gx;
	}

	public double getGy() {
		return gy;
	}

	public double getNu() {
		return nu;
	}

	private NSDiscretization setDx(double dx) {
		this.dx = dx;
		this.dxdy = this.dx * this.dy;
		this.dx2 = dx * dx;
		this.dx2dy2 = this.dx2 * this.dy2;
		return this;
	}

	private NSDiscretization setDy(double dy) {
		this.dy = dy;
		this.dxdy = this.dx * this.dy;
		this.dy2 = dy * dy;
		this.dx2dy2 = this.dx2 * this.dy2;
		return this;
	}

	public NSDiscretization setG(double gx, double gy) {
		this.gx = gx;
		this.gy = gy;
		return this;
	}

	/**
	 * Dywergencja Pola wektorowego prędkości w punkcie (i,j)
	 */
	public double D(int i, int j, Cell[][] cells) {
		double d = (cells[i][j].u - cells[i - 1][j].u) / dx;
		d += (cells[i][j].v - cells[i][j - 1].v) / dy;
		return d;
	}

	/**
	 * Wzór na składową prędkości w kierunku wersora (1,0) o współrzędnych
	 * (i+1/2,j) w następnym kroku czasowym
	 */
	public double u(int i, int j, Cell[][] cells) {
		double du = -du2dx(i, j, cells) - duvdy(i, j, cells) - dfidx(i, j, cells);
		du += nu * (d2udx2(i, j, cells) + d2udy2(i, j, cells));
		du += gx;
		return cells[i][j].u + du * dt;

	}

	/**
	 * Wzór na składową prędkości w kierunku wersora (0,1) o współrzędnych
	 * (i,j+1/2) w następnym kroku czasowym
	 */
	public double v(int i, int j, Cell[][] cells) {
		double dv = -dv2dy(i, j, cells) - duvdx(i, j, cells) - dfidy(i, j, cells);
		dv += nu * (d2vdx2(i, j, cells) + d2vdy2(i, j, cells));
		dv += gy;
		return cells[i][j].v + dv * dt;
	}

	private final double pow2(double arg) {
		return Math.pow(arg, 2);
	}

	private final double dv2dy(int i, int j, Cell[][] cells) {
		double dv2 = pow2((cells[i][j].v + cells[i][j + 1].v) / 2);
		dv2 -= pow2((cells[i][j].v + cells[i][j - 1].v) / 2);
		return dv2 / dy;
	}

	private final double du2dx(int i, int j, Cell[][] cells) {
		double du2 = pow2((cells[i][j].u + cells[i + 1][j].u) / 2);
		du2 -= pow2((cells[i - 1][j].u + cells[i][j].u) / 2);
		return du2 / dx;
	}

	private final double duvdx(int i, int j, Cell[][] cells) {
		double duv = ((cells[i][j].u + cells[i][j + 1].u) / 2) * ((cells[i][j].v + cells[i + 1][j].v) / 2);
		duv -= ((cells[i - 1][j].u + cells[i - 1][j + 1].u) / 2) * ((cells[i][j].v + cells[i - 1][j].v) / 2);
		return duv / dx;

	}

	private final double duvdy(int i, int j, Cell[][] cells) {
		double duv = ((cells[i][j].u + cells[i][j + 1].u) / 2) * ((cells[i][j].v + cells[i + 1][j].v) / 2);
		duv -= ((cells[i][j].u + cells[i][j - 1].u) / 2) * (cells[i][j - 1].v + cells[i + 1][j - 1].v / 2);
		return duv / dy;
	}

	private final double dfidy(int i, int j, Cell[][] cells) {
		double dfi = cells[i][j + 1].fi - cells[i][j].fi;
		return dfi / dy;
	}

	private final double dfidx(int i, int j, Cell[][] cells) {
		double dfi = cells[i + 1][j].fi - cells[i][j].fi;
		return dfi / dx;
	}

	private final double d2vdx2(int i, int j, Cell[][] cells) {
		double d2v = cells[i - 1][j].v - 2 * cells[i][j].v + cells[i + 1][j].v;
		return d2v / dx2;
	}

	private final double d2udx2(int i, int j, Cell[][] cells) {
		double d2u = cells[i - 1][j].u - 2 * cells[i][j].u + cells[i + 1][j].u;
		return d2u / dx2;
	}

	private final double d2vdy2(int i, int j, Cell[][] cells) {
		double d2v = cells[i][j - 1].v - 2 * cells[i][j].v + cells[i][j + 1].v;
		return d2v / dy2;
	}

	private final double d2udy2(int i, int j, Cell[][] cells) {
		double d2u = cells[i][j - 1].u - 2 * cells[i][j].u + cells[i][j + 1].u;
		return d2u / dy2;
	}

	public double fi(int i, int j, Cell[][] cells) {
		double fi = (cells[i - 1][j].fi + cells[i + 1][j].fi) * dy2;
		fi += (cells[i][j + 1].fi + cells[i][j - 1].fi) * dx2;
		fi += cells[i][i].z * dx2dy2;
		return fi / (2 * (dx2 + dy2));
	}

	/**
	 * Wyliczanie rozkładu źródeł ciśnienia
	 */
	public double z(int i, int j, Cell[][] cells) {
		return z1(i, j, cells) + z2(i, j, cells) + z3(i, j, cells) + z4(i, j, cells);
	}

	private double z1(int i, int j, Cell[][] cells) {

		double z1 = 0;

		if (cells[i + 1][j].flag != Cell.Flag.BOUNDARY) {

			z1 += (pow2((cells[i][j].u + cells[i + 1][j].u) / 2) - pow2((cells[i - 1][j].u + cells[i][j].u) / 2)) / dx2;

		}

		if (cells[i - 1][j].flag != Cell.Flag.BOUNDARY) {

			z1 += (pow2((cells[i - 2][j].u + cells[i - 1][j].u) / 2) - pow2((cells[i - 1][j].u + cells[i][j].u) / 2))
					/ dx2;

		}

		if (cells[i][j + 1].flag != Cell.Flag.BOUNDARY) {

			z1 += (pow2((cells[i][j].v + cells[i][j + 1].v) / 2) - pow2((cells[i][j - 1].v + cells[i][j].v) / 2)) / dy2;

		}

		// TODO W ksiazce jest cells[i-1][j] Traktuje to jako blad. Do
		// sprawdzenia pozniej
		if (cells[i][j - 1].flag != Cell.Flag.BOUNDARY) {
			z1 += (pow2((cells[i][j - 1].v + cells[i][j - 2].v) / 2) - pow2((cells[i][j - 1].v + cells[i][j].v) / 2))
					/ dy2;
		}

		return z1;
	}

	private double z2(int i, int j, Cell[][] cells) {
		double z2 = (cells[i][j].u + cells[i][j + 1].u) * (cells[i][j].v + cells[i + 1][j].v);
		z2 += (cells[i - 1][j].u + cells[i - 1][j - 1].u) * (cells[i][j - 1].v + cells[i - 1][j - 1].v);
		z2 -= (cells[i][j].u + cells[i][j - 1].u) * (cells[i][j - 1].v + cells[i + 1][j - 1].v);
		z2 -= (cells[i - 1][j].u + cells[i - 1][j + 1].u) * (cells[i][j].v + cells[i - 1][j].v);
		return z2 / (2 * dxdy);
	}

	private double z3(int i, int j, Cell[][] cells) {
		return -cells[i][j].div / dt;
	}

	private double z4(int i, int j, Cell[][] cells) {

		double z4 = 0;

		if (cells[i + 1][j].flag != Cell.Flag.BOUNDARY) {

			z4 += (cells[i + 1][j].div - cells[i][j].div) / dx2;

		}

		if (cells[i - 1][j].flag != Cell.Flag.BOUNDARY) {

			z4 += (cells[i - 1][j].div - cells[i][j].div) / dx2;

		}

		if (cells[i][j + 1].flag != Cell.Flag.BOUNDARY) {

			z4 += (cells[i][j + 1].div - cells[i][j].div) / dy2;

		}

		// TODO W ksiazce jest cells[i-1][j] Traktuje to jako blad. Do
		// sprawdzenia pozniej
		if (cells[i][j - 1].flag != Cell.Flag.BOUNDARY) {

			z4 += (cells[i][j - 1].div - cells[i][j].div) / dy2;

		}

		return -z4 * nu;
	}

}
